clear
m=1;
hb=1;
L=10;
dx=0.1;
x=-L:dx:L;
dt=0.001;
t=0:dt:1;
psi=zeros(length(x),length(t));%时间位置矩阵
x0=0;
Dx=1;
psi0=exp(-x'.^2/Dx^2);%初始值
psi(:,1)=psi0;%为第一列赋予初值
%引入二阶导算符的系数矩阵
A=-2*eye(length(x))+diag(ones(1,length(x)-1),1)+diag(ones(1,length(x)-1),1)+diag(ones(1,length(x)-1),-1);
%设置v势能
v=x;
V=diag(v);
H=hb^2/2/m/dx^2*A+V;
for n=1:length(t)-1
   psi(:,n+1)=psi(:,n)+1/1i/hb*H*psi(:,n)*dt;
   %两个边界条件
   psi(1:n+1)=0;
   psi(end:n+1)=0;
   plot(x,real(psi(:,n+1)),x,imag(psi(:,n+1)));
   axis([x(1) x(end) -1 1])
   getframe;
end

